Machine learning molecular dynamics reveals the structural origin of the first sharp diffraction peak in high-density silica glasses

The first sharp diffraction peak (FSDP) in the total structure factor has long been regarded as a characteristic feature of medium-range order (MRO) in amorphous materials with a polyhedron network, and its underlying structural origin is a subject of ongoing debate. In this study, we utilized machine learning molecular dynamics (MLMD) simulations to explore the origin of FSDP in two typical high-density silica glasses: silica glass under pressure and permanently densified glass. Our MLMD simulations accurately reproduce the structural properties of high-density silica glasses observed in experiments, including changes in the FSDP intensity depending on the compression temperature. By analyzing the simulated silica glass structures, we uncover the structural origin responsible for the changes in the MRO at high density in terms of the periodicity between the ring centers and the shape of the rings. The reduction or enhancement of MRO in the high-density silica glasses can be attributed to how the rings deform under compression.


Machine learning potential of silica
In this appendix, we show the details of the present machine-learning potential (MLP) and the validation of MLP for crystal structures of silica.

Details of the MLP in the present work
The reference data of MLP of silica was generated by DFT calculation using Vienna Ab initio Simulation Package (VASP).The quartz, cristobalite, tridymite, stishovite, amorphous, and liquid structures were used for constructing the DFT reference data.First, we generated various configurations for the above structures using MD simulations with Tersoff potential 1 .The energies and forces for the configurations were then reevaluated by DFT calculation and used as the DFT reference data.We employed the strongly constrained and appropriately normed meta-GGA exchange-correlation functional for the DFT calculation with an energy cutoff of 500 eV and k-spacing 0.25 Å−1 .We also conducted DFT-NPT simulations for the quartz, cristobalite, tridymite, stishovite, and liquid structures at 300, 500, 1000, 3000, and 5000 K.The DFT-NPT simulations for the crystal structures at 5000 K were performed until the solids melted completely.We picked up the structures from the MD trajectories at 25 fs intervals, and the total number of the DFT reference data for making the MLP was 21338 structures.90% of the DFT reference data was assigned to training data and the remaining 10% to test data, respectively.We used the n2p2 code 2,3 for training MLP based on the Behler-Parrinello type neural network 4,5 .We adopted the following symmetry functions as the descriptors of the distances and the angles of atoms, respectively, i.e., with the cutoff function 2 where R ij is the distance between the i-th and j-th atoms, and θ ijk is the angle between i-j and i-k atom bonds .The cutoff radius R c for G (R) i and G (A) i were taken as 8.0 Å and 6.5 Å, respectively.The other parameters of symmetry functions ( η (R) , η (A) , R s , λ, and ξ ) were selected by CUR decomposition 6 .We used two hidden layers with softplus activation functions with 20 nodes.The root mean square errors (RMSE) of energy and force for the training and test data were summarized in table S1.In the main text, we have shown that the present MLP of silica well reproduces the structural properties of silica glass.Here, we show the performance of the MLP for the crystal phase of silica.Table S2 shows the lattice constants of various crystal phases obtained by structural optimization using DFT and MLP.The lattice constants computed by DFT based on SCAN are in quite good agreement with the experimental data.The present MLP is able to reproduce the DFT results very well, with the maximum error for DFT data being only 1.43 for α-cristobalite.For the validation of our present MLP, we also compute the finite temperature properties of quartz with 1125 atoms.Firstly, we evaluate vibrational density of state (VDOS) at finite temperature using the Fourier transform of the velocity autocorrelation function as where ω, T , N , k B , m j , and v j are the frequency, temperature, the number of atoms, the Boltzmann constant, the mass of the j-th atom, and the velocity of the j-th atom, respectively.To obtain the velocities of atoms at each time step from MD trajectory, we conducted MLMD-NVE simulation at 300 K with 0.5 fs time step and 200 ps simulation period.The peak positions of the VDOS obtained by MLMD agree well with the experimental data 13 as shown in Fig. S1(a).Next, we investigate the phase transition from α-to β-quartz using MLMD-NPT simulation.Fig. S1(b) shows temperature dependence of the volume of quartz per unit cell.The quartz in α-phase shows positive thermal expansion, whereas that in βphase reveals negative thermal expansion.Although the curve of volume change obtained by MLMD slightly shift to a lower temperature than the experimental curve, the values are quite close to experimental data.From the above validations, we confirmed that the present MLP enables us to conduct accurate calculation not only for the amorphous phase but also for the crystal phase of silica.
2 Additional data of simulated high-density silica glass Fig. S2(a) shows the Faber-Ziman total structure factors of the SGUP obtained by the MLMD and high-pressure XRD experiment 16 .Although some discrepancies exist with regard to the height of peaks, the MLMD simulation is able to reasonably reproduces the high-pressure XRD experimental data, including a small peak around 2.9 Å, which can not be reproduced by some classical force field 18,19 (refer to supplementary information of reference 16 ) Fig. S2(b) also shows the total structure factors of the DSG obtained by the MLMD, and experiments 17 .
Here, we compared experimental data with calculated data of the DSG3500, which has a density close to the experimental one.The results of the MLMD simulation agree well with the experimental data of the DSG created by hot compression.Fig. S3 shows the partial structure factors S αβ (Q) of the OSG, SGUP, and DSG.Overall, the first peaks of the partial structure factors in the SGUP decrease and broaden with a peak shift towards a higher scatter vector.In the case of the DSG, the first peak in the Si-Si partial structure factor significantly developed, which primarily contributes to the enhancement of the FSDP in the total structure factor.Although both PPs in the S FZ (Q) for the SGUP and DSG become sharper than that of the OSG (see the upper panel of Fig5.(a)), the contribution of each S αβ (Q) to the S FZ (Q) differs in the SGUP and DSG.As density increases, the second negative peak intensity in the S SiO (Q) behaves oppositely: the second negative peak of the SGUP increases, while that of the DSG decreases.The small peak of the XRD S FZ (Q) around 2.9 Å for the SGUP, which is absent in the DSG (see Fig. S2), can be ascribed to the increase of the second negative peak in the S SiO (Q).shows the scatter plot of the intensity of the L 1 ∩ L 2 band peak and the FSDPs height, revealing a correlation between the two.The DSGs compressed at higher temperatures are accompanied with large structural relaxation, which make more aligned cage boundary surfaces formed by rings with small distortion contributing to the enhancement of the peak intensity of the FSDP.

Additional information of the ring approximated as a cuboid
In the main text, we have presented the change in the aspect ratio l 2 /l 3 , which characterizes the deformation of rings within a ring pseudo-plane.In this supplementary material, we show the change of the ratio of thickness l 1 and average length l ⊥ = (l 2 + l 3 )/2 of the ring pseudo-plane for different ring sizes in Fig. 6(a).The ratio l 1 /l ⊥ approaches zero when the ring pseudo-plane is almost flat.The difference in l 1 /l ⊥ of the DSG for small rings compared to that of the OSG is smaller than that of the SGUP.Conversely, the difference in l 1 /l ⊥ of the DSG for larger rings tends to be larger than that of the SGUP.The ring size-dependent changes in l 1 /l ⊥ of the DSG are the same as those of l 2 /l 3 , meaning that larger rings in the DSG are more likely to deform to avoid a substantial reduction in the Si-O-Si angle, as explained in the main text.Additionally, the distribution of l 1 and l 3 calculated for each ring size is shown in Fig. 6(b) and (c).The ring size-dependent of the l 1 /l ⊥ of the DSG is the same as those of l 2 /l 3 , , meaning that larger rings in the DSG are more likely to deform to avoid a substantial reduction in the Si-O-Si angle, as explained in the main text.We also show the distribution of l 1 and l 3 calculated for each ring size in Fig. 6(b) and (c) as additional information.

Fig. S1 :
Fig. S1: (a) VDOS of α-quartz at 300 K obtained by MLMD.Blue, red, and black solid lines show the silicon, oxygen, and total VDOS, respectively.Black dashed line denotes the experimental data of inelastic X-ray scattering 13 .(b) Temperature dependence of the volume of quartz obtained by MLMD, MD with BKS potential 14 , and experiment (Exp) 15 .

2. 1 1 Fig
Fig. S2: (a) Faber-Ziman total structure factors of XRD of the SGUP.Solid lines (SGUP) are the results obtained by MLMD simulation and the dashed lines (EXP1 and EXP2) are experimental data taken from supplementary information of reference 16 .(b) Faber-Ziman total structure factors of ND and XRD of the DSG.Solid lines (DSG3500) is the result computed by MLMD simulation and the dashed lines (EXP3 and EXP4) are experimental results 17 .The values in (•) in the legends of (b) are the density of silica glass [g/cm 3 ].

2. 2 Fig. S3 :
Fig. S3: Partial structure factors for the O-O, Si-O, and Si-Si pair.The black, blue, and red lines are the structure factors for the OSG, the SGUP, and the DSG, respectively.The values in (•) in the legends are the density of silica glass [g/cm 3 ].

2. 3
Fig. S4: (a) and (b) Ring center differential distribution function GRC(r) for OSG, SGUP, and DSG3500.(c) and (d) Ring center structure factor SRC(Q) for OSG, SGUP, and DSG3500.The values in (•) in the legends are the density of silica glass [g/cm 3 ].

FigureFig
FigureS4shows systematic changes in the ring center differential distribution functions G RC (r) and structure factor S RC (Q) of the SGUP and DSG3500.With increasing density, the G RC (r) peaks of the SGUP in the intermediate range move towards shorter distance, and the peaks of S RC (Q) shift to a higher scattering vector Q.The shapes of the S RC (Q) are almost identical to that of the OSG, indicating that the origin of the MRO in the SGUP is the same as that in the OSG, as discussed in the main text.In contrast, the DSGs show a new peak at around 12 Å, in addition to the shifts of the medium-distance peaks in G RC (r).Furthermore, the S RC (Q) peaks become sharper as shown in Fig.S4(d), indicating the development of quasi-periodicity in the DSGs.

Fig
Fig. S6: (a) Averaged ratio of l1/l ⊥ for each ring size.(b) and (c) show the distribution of l1 and l3 calculated for each ring size, respectively.The values in (•) in the legends are the density of silica glass [g/cm 3 ].

Table S1 :
RMSE of MLP for the training and test data.

Table S2 :
The lattice constants of various crystal structures of silica obtained by DFT, MLP, and experiment (EXP).